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Abstract 

We have studied the presence of plateaus on the low-temperature magnetization of an antiferro¬ 
magnetic spin-1 chain, as an external uniform magnetic held is varied. A crystal-held interaction is 
present in the model and the exchange constants follow a random quenched (Bernoulli or Gaussian) 
distribution. Using a transfer-matrix technique we calculate the largest Lyapunov exponent and, 
from it, the magnetization at low temperatures as a function of the magnetic held, for different 
values of the crystal-held and of the width of the distributions. For the Bernoulli distribution, the 
number of plateaus increases, with respect to the uniform case (F. Litaiff, J. R. de Sousa, and N. 
S. Branco, Sol. St. Comm. 147, 494 (2008)) and their presence can be linked to different ground 
states, when the magnetic held is varied. For the Gaussian distributions, the uniform scenario is 
maintained, for small widths, but the plateaus structure disappears, as the width increases. 
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I. INTRODUCTION 


The physics of one-dimensional or quasi one-dimensional, spin-S' (S' ^ 1), chains have 


fl 


attracted a considerable amount of attention since Haldane’s prediction [fj that the one¬ 
dimensional Heisenberg antiferromagnetic model should have an energy gap between the 
singlet ground state and the first excited triplet states in the case of an integer spin quantum 
number, while the energy levels are gapless in the case of a half-integer spin quantum number. 

One fascinating characteristic of some low-dimensional quantum systems is that they 
show magnetic plateaus, i.e., quantization of the magnetization at low temperatures near 
the ground state. A general condition for the quantization of the magnetization in low¬ 
dimensional magnetic systems is derived from the Lieb-Schultz-Mattis theorem. One 
important issue on this matter is the effect of a single-ion anisotropy term on the critical 


and multicritical phenomena, which 


such as, mean-field approximation jsj, transfer-matrix, finite-size-scaling and Monte Carlo 
simulation jd-T], Bethe lattice treatment [^, and effective held theory |l^, Q. The 
mechanism for the appearance of these magnetization plateaus in quasi one-dimensional spin 
chains are expected to be, among others, dimerization, frustration, single-ion anisotropy, 
quantum effects, and periodic held. However, plateau structures may also be present in 
classical models at zero temperature. These structures, as we will see, are connected to 
diherent ground-state conhgurations. 

Experimentally, the magnetization plateaus were observed in high-helds measurements 


las been investigated by various theoretical methods. 


of several magnetic materials such as the quasi one-dimensiona 


Y 2 BaNiO^ ll3|, Ni{C2HsN2)2N02Cl0^ (abreviated NENP) 
16l |: triangular antiferromagnets CqEu 17|, CsCuCl^ 
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compounds SrCu 20 ^ |l2| . 


1^, and Cu{N 0 ^) 22 .bH 20 
and RbFe{M 004)2 19|; 


and in a quasi two-dimensional compound, with a Shastry-Sutherland lattice structure, 
SrCu 2 {BO ^)2 |20|. The magnetic susceptibility of this compound shows a couple of charac¬ 
teristics of a spin gap system: a Curie-Weiss behavior at high temperatures, with a round 
maximum at 20iC, followed by a sudden decrease and zero spin susceptibility at T = 0 


2l| . Another piece of evidence for the extremely localized nature of triplet excitations is 


obtained when one measures the magnetization curve at low temperatures. Onizuka et ah 
^ observed a predicted plateau at 1/3 of the total magnetization around 50 T, in which 
the magnetic superstructure is characterized by a novel stripe order of triplets. He observed 


3 































that, irrespective of the direction of magnetic held, the magnetization of SrCu 2 {BOz )2 is 
quantized at 1/8, 1/4, and 1/3 of the fully saturated Cu moment. 

Our goal in this work is to model the presence and conhguration of plateaus in a one¬ 
dimensional system with quenched disorder. The presence of impurities is ubiquitous in 
nature and the understanding of their effect is, therefore, of great theoretical and experi¬ 
mental relevance. Here we focus on random and frozen impurities, as opposed to determin¬ 
istic inhomogeneities (like aperiodic modulation, for example); the latter will be studied in 
a forthcoming work. 

This work is organized as follows. In Section HIl we present the model and formalism. 
In Section IIIII the numerical results and discussions are given. Finally, the last section is 
devoted to conclusions. 


II. MODEL AND FORMALISM 


In general, most of the theoretical studies on low-dimensional spin-S' systems have been 
based on Heisenberg Hamiltonians, which presents non-trivial quantum effects. However, 
several studies have shown that it is possible to use classical spin systems to obtain magnetic 
plateaus and study the dependence of these on physical characteristics |23l-l28|. as discussed 


in the Introduction. Consider, for example, the Hamiltonian of a one-dimensional spin-1 
classical system, which uses the Ising type variable instead of quantum spin operators: 


« = J ^ S'SJ., - J/ ^ Sf + A )^ 

i i i 


( 1 ) 


where J stands for the antiferromagnetic (J > 0) exchange coupling, A is the single-ion 
anisotropy (or the crystal-held term), H is an external uniform magnetic held, and Sf = 


± 1 , 0 . 


In a previous work 2^, it has been shown that this model allows for magnetization 

ate thermodynamic quantities for 


plateaus at zero temperature. The method used to calcu 
this Hamiltonian [2^ was the transfer-matrix technique 


29 


30| . which, besides being very 


simple, allows for exact results when applied to one-dimensional classical systems. Within 


this formalism, it can be shown 
the transfer matrix A4, by: 


3l[ | that the partition function Z is given, as a function of 
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Z = Tr>l^ = Af+ A^ + Af+ ... ~Af, (2) 

ioT N ^ oo, where N is the number of spins in a one-dimensional lattice or the number of 


rows in a two-dimensional one and Ai > A 2 > A 3 > ... are the eigenvalues of Ai. 

Therefore, the free energy per spin reads 

fiH,T) = -ksT lim -fcsThiAi, (3) 

V-foo iV 

where /c^ is the Boltzmann constant and T is the temperature. 

For the model dehned by the Hamiltonian in Eq. ([T]), the transfer matrix reads 


where x = exp(—JF), y 
^ = l/kBT. 


M 


^ xyz yz yz/x'^ 
1 1 1 
\y/{.xz) y/z xy/z j 


exp(—d) and z = exp(L), with K 


/3J, d = PA, L 


(4) 

pH, and 


Our goal in this work is to study the effects of disorder in the plateaus’ structure. We 
introduce the disorder in the exchange constant J, which will have its value given by a 
probability distribution (see below). Therefore, the Hamiltonian now reads 


w = 5] S' + A 5](S')^ 

i i i 

such that the transfer matrix is given by 




yz 

yz/xi ^ 

-M,= 

1 

1 

1 



y/z 

Xiy/z j 


with Xi = exp(—iFj) and Ki = pji and all other dehnitions as in Eq. (jl]). 
We have used two different probability distributions, namely: 


(5) 


( 6 ) 


V{JP = V J {1 + a)) + il-p) J {I - a)), 


where p = 1/2 and d is the Dirac delta function, and 


V(Ji) 


1 


\j2'Ka‘^ 


exp 


2a2 


(7) 

( 8 ) 
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In both cases, the width of the distribution, in unities of J, is given by a. 

The expression shown in Eq. (3) is in fact a simplihcation suitable for uniform systems. 
On more general grounds, the partition function is obtained from applying the transfer 
matrix to a suitable initial vector a great number of times, such that the dominant behavior 
is obtained (this behavior is connected to the greatest eigenvalue of the matrix, for uniform 
models). Therefore, the correct expression for the calculation of the free energy per site in 
the thermodynamic limit and for a disordered model, within the transfer-matrix framework, 
reads: 


ksT 


—Ai = — lim In 

N^OO 



(9) 


where now the transfer matrix Adj is a stochas 
a given probability distribution (see Eq. |6]) [32 


ic quantity and its elements are taken from 


3^, is the size of the strip, ho is a suitable 


vector to start the iteration with, and the disordered transfer matrix is applied N times to 
ho- This starting vector was varied and the results do not depend on it, for iV S> 1. In the 
previous expression, Ai is the greatest Lyapunov exponent 3^. It is readily seen that Ai is 


the greatest eigenvalue of the transfer matrix for uniform models. 

Eq. (I9j) holds for a given realization of the disorder; since we are dealing with a quenched 
model, the free energy has to be averaged for different realizations of the disorder. This 
quenched average is represented by [...] and, then, we have 


'ave 


ksT 

such that the magnetization is given by 


-[Ai], 

Q jave 

OH ’ 


( 10 ) 


( 11 ) 


where the derivative is numerically obtained. 

Some technical points are worth stressing here. The hrst relevant point is the size of the 
chain, N, which is supposed to be big enough such that, after successive applications of the 
transfer matrix, a given accuracy is obtained for the free energy. For practical purposes, we 
have iterated Eq. ([9]) until the values for the free energy do not differ by more than 0.1%; 
roughly, this precision is achieved for N > 10^. For this value of N^ the number of disorder 
realizations, Nd, which lead to an error of 1% or less in the magnetization is Nd ~ 10. The 
error in is obtained as one standard deviation among the Nd values. For the numerical 
derivatives involved in Eq. flTT]) . a value of 5H (in units of J) of order 10“^ was used. 


6 













Another important point is that onr calculation is made at low temperatures. Therefore, 
the interaction constant K is big and the iterated vectors may have big components. To avoid 
overflow, we normalize the vectors after each iteration of the transfer matrix and use their 
norms to calculate the free energy, using Eq. ([2]). Note also that, for finite temperatures, the 
free energy is an analytic function of the field. This analyticity precludes mathematically 
flat regions, with sharp transitions. However, on the scale of measurements or on numerical 
calculations, such plateau regions can and do occur. Moreover, an extrapolation procedure 
to the zero-temperature limit clearly shows that the usual plateau structure, with sharp 
transitions between steps, is obtained. 


III. RESULTS AND DISCUSSION 


First, we show our results for the uniform case, in order to test our procedure and code. 
The results are presented in the inset of Fig. HI where the magnetization m is depicted as a 
function of the reduced magnetic held, H/J, for A/J = 1: they agree with those obtained 


in Refs. 


23 and 


26l . In fact, since our procedure is less time- and memory-consuming, we 


are able to go to temperatures closer to zero and there is less rounding in our plateaus than 


for the ones in Ref. 


23 


For 0 < A/J < 1, the hrst plateau (see inset of Fig. [T]), at m = 0, is linked to the 
conhguration (h—), the second one, at m = 1/2, rehects the conhguration (-1-0), and the 
last one is connected to the conhguration (-I--I-). Here, -1-, — , and 0 stand for the values 
-|-1, —1, and 0 for the spins, respectively. The actual conhguration of the whole chain is 
obtained by repeating the two-spin structure inside parentheses. The hrst step happens at 
H / J = 2 — A/ J and the second one takes place atiJ/J = 2-|-A/J. On the other hand, for 
A/J > 1, the hrst plateau, with m = 0, is connected to the conhguration (00), while the 
other two plateaus correspond to the same conhgurations as for 0 < A/J < 1. The step at 
the left-hand side now happens at H/ J = A/ J and the step at the right-hand side is placed 
in the same position as for 0 < A/J < 1. We have also calculated the magnetization vs 
magnetic held behavior for A/J < 0: only two plateaus are present, since negative values 


of A do not favor the zero state 


26|. 


We now present our results for the disordered model given by the Bernoulli distribution 
(see Eq. ([7])). In Fig. [H where the plateaus structure is shown for three diherent values 
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of the width a and for A/J = 1 (the qualitative results do not differ for other values of 
A/J > 0; for A/J < 0, as stated above, there are only two plateaus for the uniform model 
but the introduction of disorder brings the same overall structural changes as for the three- 
plateau conhguration), for a temperature T = 0.01 (henceforth, we refer to the temperature 
in unities of J/ks, i.e., ksT/J). We notice that the number of plateaus increases with the 
disorder but we still can connect them to different ground-state conhgurations. 

In order to make this connection, we depict, in Fig. [2l a portion of the chain, with 
weak and strong interactions represented by thin and thick lines, respectively. Due to the 
random nature of the interactions spins may be connected to their two hrst-neighbors by 
weak interactions only, by strong ones only or by one of each. In what follows, we show 
some examples that link the steps in Fig [T] with ground-state conhgurations. 

Let us imagine that the magnetic held is strong enough that all spins are in the -|-1 state; 
this corresponds to m = 1 (see Fig. [1]). When the held is decreased, spins connected to their 
two hrst neighbors by strong links will change hrst, to the 0 state. But not all spins with this 
neighborhood will turn, since these changes increase the contribution of the magnetic-held 
term to the energy (although decreasing the corresponding contribution of the exchange and 
crystal-held terms to the Hamiltonian). It is energetically favorable to turn every other spin 
in a cluster of contiguous strong interactions. In Fig. [3] (a), spins represented by S 2 and 5*4 
(or by S 3 and S^) will turn hrst, when H is decreased, while in Fig. |3] (fe) spins S 2 and S '4 
will go to the 0 state when the magnetic held is decreased. In any case, the quoted spins 
will turn when the magnetic held is decreased and crosses the critical value He'. 


He/J = 2 (1 + a) + A/J, 


( 12 ) 


which correspond to the three rightmost steps in Fig. [T] (recall that A/J = 1). To calculate 


the drop in the magnetization on these steps, one has to work out 


islands of bonds, which is a one-dimensional percolation problem 


)he statistics of random 
. The probability of 


a cluster of s strong bonds to be generated is p^(l — p)^, while the number of spins that 
change state is Ug = s/2, for even s, and Ug = (s — l)/2, for odd s. Therefore, the fraction 
of spins that change state, is: 


= 




s=2,s even 


(13) 


s=3,s odd 


where the hrst sum is over even values of s and the second one over odd values of s. 




For p = 1/2, we obtain ps = 1/6 = 0.166. Therefore, the magnetization drops to 

5/6 = 0.833..., independent of a. 

At H/J = 3 another step is present, now independent of a, and the drop in m is again 
1/6. The spins that change state are those in clusters such that the spins interact with their 
nearest neighbors via one strong and one weak bond, like the ones depicted in Fig. 01 The 
spins which change state are, for example, S 2 and 5*4 (or S '3 and S^) in Fig. |4] (a) and 5 * 2 , 5 * 4 , 
and Fe in Fig. 01(6). The calculation is the same as for the preceding case and, therefore, the 
drop on the magnetization is also 1/6. However, the energy balance is a little bit different 
and the step takes place at 

HJJ = 2 + A/J. (14) 

The next step is connected to spins on islands of contiguous weak bonds and the pro¬ 
cess is the same as the one described above for islands of strong contiguous bonds. The 
magnetization decreases by 1/6 and the steps are located at: 

HJJ = 2{1-a) + A/J (15) 

Equivalent steps are obtained when spins on the 0 state turn to the —1 state. These steps 
are placed on smaller values of H/ J. 

Therefore, the presence of random disorder described by a Bernoulli distribution increases 
the total number of steps, when compared to the uniform model, for the same valeu of A/ J. 
However, the steps are still connected to the energy balance between different conhgurations 
of ground states (see Fig. [1]). 

We now turn to a Gaussian distribution. We calculate the magnetization as a function 
of A/J, the magnetic held H and the variance a of the Gaussian probability distribution 
given by Eq. (0). The numerical results are obtained as follows: in Figs. [5] and El the 
magnetization m is plotted at T = 0.01 as a function oi H/J and for four values of the 
anisotropy parameter (namely, A/J = 0.5, 1.0, 1.5, and 2.0), with a = 0.1 and a = 0.2 
respectively. 

As expected (see Fig. E]), the results for low a are qualitatively equivalent to those for the 
uniform model and three plateaus are present for A/J > 0. This equivalence comes from 
the fact that, for small values of a, the Gaussian distribution is concentrated around the 
value Ji = J (see Eq. ([ 8 ])). This mimics the uniform model, where all exchange interactions 
are the same. For the Bernoulli distribution, on the other hand, no matter the (hnite) 
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width of the distribution, the two peaks have the same height. Therefore, two different 
exchange parameters will be involved in this model. This behavior is depicted in Fig. [T] 
as long as (T 7 ^ 0, there will be more steps than for the uniform model. But we can see 
that some steps will disappear when a goes to 0, regaining the uniform model behavior. 
Note, however, that the width of the plateau diminishes when disorder is introduced, with 
the corresponding rounding of the steps. This is consistent with the interpretation that 
the steps are connected to transitions between different ground states. The continuous 
probability distribution leads to a variety of strengths for the exchange interaction, such 
that different spins change state at different values of H/J. Note that the steps are still 
centerd at the values for the corresponding steps for the uniform model (as discussed in the 
hrst paragraph of Sec. HID; the effect of the distribution is only to smooth some of these 
steps. Also, the left-hand side steps for A/J = 1.5 and 2.0 are not rounded; the reason 
is that the energy balance for the transition between states (00) and (-1-0) (see again hrst 
paragraph of Sec. IIIIll does not involve the exchange constant. 

The picture discussed in the last paragraph is supported by the results depicted in Fig. |6], 
for (T = 0.2. A wider distribution for J allows for a more pronounced rounding of some steps 
(the left-hand steps for A/J = 1.5 and 2.0 remain as for the uniform model, as discussed). 
Another consequence of this picture is that, for A/J = 0.5, the intermediate plateau is 
not present anymore (compare with Fig. [5l where this plateau, for A/J = 0.5, still exists, 
although it is narrower than for A/J > 1). 

In Figs. [7] and [HI the magnetization m at T = 0.01 is plotted as a function of H/J for 
several values of the variance a (0.1, 0.2, and 0.3) with A/J = 1.0 and 2.0, respectively. 
These hgures show that the steps are rounded when the width of the Gaussian distribution 
is increased. 

Our results indicate that the Gaussian distribution only rounds the steps for small values 
of its width. For large ones, the continuous distribution tends to eliminate the steps and 
plateaus. This is consistent with the observation that these plateaus are consequence of 
transitions between different spin conhgurations at the ground state. 

The presence of magnetic plateaus in the present system conhrms that this structure may 
□e obtained in classical models. Similar behavior has been observed in theoretical studies 


23j, 


37| and in experimental results 15j. Note that the ground state phase diagrams are 


symmetric with respect to H/J = 0 line, i.e. invariant under the transformation H —)■ 
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—H and S —?> —S. As a final note, we would like to mention that the effects of disorder 
depend on the model and may be completely different from the ones observed in this work. 
When applied to a quantum model, for example, disorder following a Gaussian distribution 
completely eliminates the presence of plateaus 


IV. CONCLUSIONS 


In this work, we studied the low temperature magnetization plateaus of a classical spin 
system with discrete symmetry. By using a transfer-matrix technique, we calculate the 
magnetization behavior at low temperatures as a function of the anisotropy parameter A/J, 
magnetic field H, and variance a of the probability distributions. It is found that positive 
single-ion anisotropies have a significant effect on the low-temperature magnetic properties 
in the system, as the external held is varied. The number of plateaus for the Bernoulli 
distribution increases, when compared to the uniform model, and their presence is linked to 
transitions between different ground states. For the Gaussian distribution and small width a, 
only one intermediate plateau (at m = 1/2) is found. The steps in the magnetization-versus- 
magnetic held graphs are rounded as the width of the Gaussian distribution is increased and 
the overall picture is consistent with the presence of steps due to transitions between diherent 
ground states. Note, however, that the ehects of the introduction of a Gaussian distribution 
depends strongly on the underlying model (see, for example. Ref. [3^. Our results are 
compatible with previous experimental and theoretical studies. Work is under way on spin- 


alternating chain models, which are expected to show many interesting phenomena 


39]. 
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0 1 2 3 4 5 

H/J 


FIG. 1: Low-temperature magnetization plateaus for the disordered model defined by the Hamil¬ 
tonian in Eq. (I5|), with A/J = 1. The exchange constant is obtained from a Bernoulli distribution 
(Eq. ([7|)) with width a. In the inset, we show the plateau structure for the uniform model, for 
comparison. 


• • • a 



EIG. 2: In this figure we represent a possible portion of the one-dimensional lattice. Weak interac¬ 
tions are represented by thin lines (like the one between spins S 2 and S 3 , for example) and strong 
interactions are represented by thick lines (like the one between spins S 3 and S 4 , for example). 
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S3 S, 


(a) 


S5 s, 


-o 

S, 


(b) 

FIG. 3: Here we depict clusters of spins connected by strong interactions with an even number 
of spins (a) and an odd number of spins ( 6 ). In the former, spins S 2 and ^4 (or S 3 and ^ 5 ) will 
change first to the 0 state when H is decreased. In the latter, spins S 2 and S '4 will turn first under 
the same conditions. 


3 

(a) 


S. 


-o 

S7 


(b) 


FIG. 4: Here we depict clusters of spins which interact with their nearest neighbors via one strong 
and one weak bond with an even number of spins (a) and with an odd number of spins ( 6 ). In the 
former, spins S 2 and S '4 (or S 3 and ^ 5 ) will change first to the 0 state when H is decreased. In the 
latter, spins S 2 , S^, and Sq will turn first under the same conditions. 
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FIG. 5: Magnetization versus magnetic field for different anisotropy parameters A/J, with variance 
a = 0.1, and at T = 0.01, for the Gaussian distribution, Eq. ([8|). 



FIG. 6: Magnetization versus magnetic field for different anisotropy parameters A/J, with variance 
a = 0.2, and at T = 0.01, for the Gaussian distribution, Eq. (I8|). 
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FIG. 7: Magnetization versus magnetic field for three different values of cr, with anisotropy param¬ 
eter A/J = 1.0, and at T = 0.01, for the Gaussian distribution, Eq. Q. 



FIG. 8: Magnetization versus magnetic field for three different values of cr, with anisotropy param¬ 
eter A/J = 2.0, and at T = 0.01, for the Gaussian distribution, Eq. ([8]). 


17 




























